program main

    use rkf45_module, only: rkf45, rk
    implicit none
    integer, parameter :: neqn = 1
    real(rk) :: t = 0.0_rk, tout = 1.0_rk, y(neqn)
    real(rk) :: relerr = epsilon(1.0_rk), abserr = epsilon(1.0_rk)
    integer :: flag = 1
    integer :: iwork(5), i
    real(rk) :: work(3 + 6*neqn), steps(14)

    
    y = 0.0_rk
    steps = [real(rk) :: 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.1, 1.2, 1.3, 1.35, 1.4, 1.45, 1.50, 1.55]
    
    print *, "   xn             yn          y(xn)=tan(xn)  y(xn)-yn"
    print '(*(F12.7, 2x))', steps(1), y(1), fx(steps(1)), y(1) - fx(steps(1))
    
    do i = 1, 13

        call rkf45(fcn, neqn, y, steps(i), steps(i + 1), relerr, abserr, flag, work, iwork)
        print '(*(F12.7, 2x))', steps(i), y(1), fx(steps(i)), y(1) - fx(steps(i))

    end do

contains

    !> Evaluates the derivative for the ODE
    pure subroutine fcn(t, y, yp)
        real(rk), intent(in) :: t
        real(rk), intent(in) :: y(:)
        real(rk), intent(out) :: yp(:)

        yp(1) = 1.0_rk + y(1)**2

    end subroutine fcn

    !> Exact solution
    pure real(rk) function fx(t)
        real(rk), intent(in) :: t

        fx = tan(t)

    end function fx

end program main
